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I. INTRODUCTION 


iiemwho 7s ,OmIumCcriIsitS. where the price of a barrel of oil 
jumped higher than 210 percent and affected every aspect of 
the world economy,directly affected the shipping transporta- 
tion industries.Whereas the price of oil was formerly | of 
modest importance, it became a prime concern.ExtenSive asso- 
ciation to the emerged problem has’ to do with the design of 
autopilots for ships, which can minimize propulsive losses 
caused by added resistance due to steering. 

One motivation for the design of an optimum steering 
controller was the work done by Nomoto and Motoyama who 
claimed that reduction in propulsion loss could amount to a 
l percent Savings in fuel consumption. 

For most commercial ships alto 2 percent saving in 
fuel costs,justifies the expense of fitting an autopilot 
which has the capability of producing this savings [Ref. 1]. 

Chapter 2’ addresses what type of computer model can be 
used to represent the ship.Several of these models are 
investigated,such as simulation from the equations of motion 
and the Nomoto third order which was developed from the 
equation of motion. 

Chapter 3 addresses the problem of selecting an adequate 
cost function to represent the added resistance due _ to 
steering. 

The problem of finding the best controller design to 
provide a minimum value of added resistance due to steering, 
for regular and irregular seas,is studied in Chapters 4 and 
5 respectively by using as a model of the ship the equations 
of motion and a function minimization subroutine. 

Chapter 6 indicates how the controller parameters can be 
adjusted in any encounter environmental condition by means 


of an adaptive control. 


dh 


Chapter 7 presents the conclusions from the experimental 


work. 


eZ 


Lives ome iio er COMPULER MODELS 


The most accurate model which can represent’ the ship/ 
steering dynamics is the model which is based upon the equa- 
tions of motion as defined by series expansion including all 
terms (both linear and nonlinear). 

Using experimentally measured hydrodynamic coefficients, 
for the SL-7 high speed containership, a computer program 
was developed in order to provide a computer simulation for 
the ship.Figure 2.1 shows the block diagram and the computer 
program can be found in Appendix A [Ref. 2,3]. 


EQUATIONS 


CONTROLLER os 
MOTION 


FUNCTION 
MINIMIZATION 
SUBROUTINE 
GEO.) 


J= (nlf 5 at 
O 





Bagptiine 2. 1 Block Diagram of Ship and Control System 
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Fa pices 22 Determination of Nomoto Third Order Model 


As a next choice, from the equations of motion we derive 
the Nomoto third order transfer function. Figure 2.2 shows 
the block diagram. 

Using the scheme of Figure 2.2 which includes the func- 
tion minimization subroutine with both yaw and Sway equa- 
tions, with the linear terms only we obtain the appropriate 
coefficients of the Nomoto third model. (Including in the 
equations of motion nonlinear termS we can see that’ the 
perturbations were small enough). 

The function minimization Subroutine used was BOXPLX, 
which was programed by R. Hilleary.The task of the already 
mentioned subroutine is to find the minimum of any function 


and is subjected to explicit constraints of the variables or 
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implicit constraints on functions of the variables.In addi- 
tion it can handle a maximum of 25 variables [Ref. 4]. 

The extracted results were very close to the results we 
got by the analytic solution and are tabulated below only 
moeelo, 23 and 32 knots. 
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III. STEERING PERFORMANCE CRITERION 


Studying the literature, one can see many approaches to 
the problem of optimizing an automatic ship’ steering 
controller for maximum reduction in fuel conSumption. Since 


added resistance is directly related to both rudder activity 


and yawing motion we can express a measure of this added 


resistance given as a performance criterion with the 
formula: 
y Ce 
1 
J= ——]( \Wo+8 Jat 
ie e (2b) 
O 
Where: 


* Delta Gopmmnade: angle 
¢ Psi (YY) =yaw angle 
° Lambda( A )=weighting factor 


This formula defines an approximate drag due to steering 
for small amplitude oscillations about a steady-state pivot 
point of the ship during yawing at the natural frequency of 
the closed-loop ship/steering system. (About 0.05 rad/sec for 
the SL-7 containership).It is also convenient for shipboard 
use because yaw error and rudder angle can be easily meas- 
ured. 

During this study the values for the weighting factors 
are taken from R. E. Reid's work for the high-speed contain- 
ership SL-7 [Ref. 5]. 
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TABLE 2 
Weighting Factors 


Ship's Speed,Knots Lambda 
16 Roe iG 
25 8.128 
32 4.2 


In Table 2 weighting factors for the operating range of 
the ship are tabulated. 
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IV. CONTROLLER DESIGN FOR REGULAR SEAS 


Our goal now, is to estimate the system's performance in 
a seaway. To accomplish the above task first we must deter- 
mine a suitable representation of the external disturbances 
on the ship by tthe sea. This can be done, Since a suffi- 
ciently accurate computer ship model and a steering perform- 
ance criterion have been defined. 

It can be postulated that a sufficiently accurate model 
of the seaway itself, is a representative modeling of forces 
and moments exerted on the ship [Ref. 6,7,8]. An explana- 
tion of what a sea comprises, and how predicted or observed 
sea states can be analyzed to determine the forces’ and 
motions of a body in a sea was studied by Michael [Ref. 9]. 

In this chapter we will uSe aS a Seaway representation 
the regular sea model in which the forces on the_~= ship 
[Ref. 10], have the form: 


Wy, Re COS( Wet + 3;) 


(4.1) 
Where: 
° Re =exciting force 
eWe =encounter frequency 
¢ Wh, =wave height 
oA =phase angle 
Values for the exciting force ( ) for different 


Emeoumeer angles and different encounter frequencies as 
well, were taken from [Ref. 3]. Table 3 shows the corre- 
spondence between wave height and sea state which were used 
in this work. 
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The controller used in the present study is’ shown in 


Appendix B and its form is repeated here in Figure 4.1 


TAREE 
Sea State Vs Wave Height 


Sea State Wave Height 
Beaufort scale in feet) 
( I ( 4.0 

2 4.5 
4 ae, 
6 10.0 
7 17.5 


K,(1+T,S) (147,38) 


(14+T7S) (1+T3S) 





Figure 4.1 Controller C 


Values for the optimal gains for the above controller 
and values for the cost J as well, are shown below in Tables 
4,5, 074 Oo OL 

These values were obtained using: 

* Sea states 1-2-4-6-7 (Beaufort scale). 

¢ Encounter angles 0°-30°-60°-90°. 

* Encounter frequencies 0.2-0.4-0.6-0.75-1.5 rad per sec. 

* Constant speed 23 knots. 


A careful analysis of the extracted results lead us to 


conclude: 
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¢ The maximum deviation of controller parameter values 
occured at 0.4 rad per sec encounter frequency for all 
used encounter angles and sea states. 

e For the same encounter frequency,for all encounter 
angles studied, the controller parameter values 
increasing smoothly as the sea state increases. 

For all encounter frequencies and for Specific 
encounter angles, we observe that the cost increases at 
higher sea states. 

e For 0.4 rad per sec encounter frequency the cost 
changes rapidly for sea state 6 and 7 as we move 
through the encounter angles. 

¢ For the same encounter angles and sea states the cost 
decreases at higher encounter frequencies. 

* For all encounter angles and encounter frequencies the 
maximum cost occured for sea state /7. 

¢* The cost had it's maximum value for encounter angle 90° 
and encounter frequency 0.4 rad per sec. 

Furthermore in order to obtain the behavior of the 
rudder and yaw motion of the ship, transient response plots 
were obtained for controller 'C' at ship's speed 23 knots, 
different sea states and encounter angles as shown in 
Figures 4.2 through 4.15 


These plots were obtained using the program of Appendix 


From these plots it is verified that for the same 
encounter frequencies and encounter angles, rudder and yaw 
motion increases in amplitude as the sea states rises. 

Finally Figures 4.16,4.17,4.18,4.19,4.20, show the 
‘results of an experiment in which we were changing one 
controller parameter keeping the rest fixed, and we observed 
that the cost does not change significantly in the vicinity 


ot the actual value. 


Zo 


Using that as a fact we can postulate that tor a 
Specific sailing mode for the controller parameter values 


high accuracy is not required. 


Comparing controller C with controller A it is obvious 


that the parameter values of controller A vary over wider 


range. 
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TABLE 6 
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Controller C For 60° Encounter Wave Angle 
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V. IRREGULAR SEAS-CONTROLLER DESIGN 


A sea State generator program which genarates added 
mass, added inertia values and in addition calculates the 
forces and momentsS,was coupled to the fortran program shown 
in Appendix D, so that the function minimization subroutine 
(BOXPLX) could be used in the presence of the irregular 
sea.The forces and moments were stored ina look up table 
which was coupled to the equation of motion. 

The optimal gains obtained for 0.75 rad/sec encounter 
frequency in the regular sea study were used as the initial 
guess in order to evaluate the optimal controller parameters 
involving the irregular sea. 

It is known that sea is never regular but actually is a 
random phenomenon where waves are continually changing in 
height, length and breadth. 

Since the sea state during this study is represented by 
irregular waves then the waves impinding on the’ ship hull 
would contain the total energy density spectrum. Thats 
dencity specrum iS composed of many frequencies and there- 
fore the response of the ship would be to an average value 
of added mass and added inertia. 

For this study where the ship's speed is 23 knots,the 
controller has the form of equation B.1, we used values for 
added mass and added inertia corresponding to 0.75 rad/sec 
because close to that frequency the energy density is 
maximum. 

The optimal controller parameters found are shown in 
Table 12,and typical system's response is shown in Figures 
5.1,5.235o eae 


These plots were obtained using the program of Appendix 
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PGArerlmanalysiS Of ithe Extracted results leads to 
conclude: 

e The maximum deviation of controller parameters values 
occured at 30° encounter angle. 

¢ For all encounter angles the maximum cost occured for 
Sea state /. 

* For specific encounter angles the higher the sea state 
the higher the cost. 
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VI. ADAPTIVE CONTROL 


A. STRUCTURE 


When the ship is moving in a seaway, the controle 
parameters are changing due to alterations in the sea state 
and encounter angle. In addition we know that using fixed 
parameters for the controller over the entire spectrum,it is 
somewhat difficult to have an appropriate response of the 
system.The adjustment of the controller parameters during 
operation in a seaway,can be achieved by means of an adap- 
tive control. 

The adaptive controller can be built with digital 
circuits and analog components as well. Analog system hard- 
ware has to be designed for each specific requirement, and 
any new requirement involves changes to components. This is 
a time consuming task.However,for simple control require- 
ments anolog systems still have a possible economic advan- 
tage over digital systems [Ref. 11]. 

On the other hand,digital systems are immediately more 
attractive when control systems are required to carry out 
more and more complex tasks.The advent of microprocessors 
and associated components has enabled low cost microcom- 
puters to be built.These no longer require special environ- 
ments and are fully compatible with shipboard use.The 
advantages of microcomputer controls over other types are 
that they are extremely reliable particularly as relatively 
fewer components are required and hence they are smaller. 
Their capability is greater than comparable analog systems 
due to their ability to carry out more complex calculations. 
A major advantage is their flexibility while using standard 


hardware, being able to be reconfigured for changes in 
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system requirement without the need to alter the _ hard- 
ware.This flexibility is achieved by the programmed software 
which is stored in the memory of the computer. 

An adaptive control scheme is indicated in Figure 6.1 
and in Figure 6.2 we show analytically the components of 


mat we call the ‘Decision Generator’. 


B. SEARCH ALGORITHM 


The parameters corresponding to various environmental 
conditions(sailing modes) are sorted in a memory.When the 
sailing modes change, delta(O ) and psi(W) vary so that the 
value of the performance index changes from its theoretical 
value. This change is detected by the threshold detector. 

If the threshold is too small the adaptive structure 
will try new conditions when it is not necessary and in the 
opposite case the adaptive structure might not adapt to new 
Sailing modes. 

Using the term threshold we want to make sure that only 
Significant changes in the environmental conditions will be 
considered, otherwise the system will look for a new model 
too often. 

When the threshold detects that a new set of parameters 
Should be used the search algorithm will look for conditions 
close to the previous one because the external conditions 
are not expected to vary rapidly. Then a new optimal value 
of J is selected and must’ be matched with the performance 
index computed for the real psi and delta. 

Matching operation must be done after some delay 
ensuring that the conditions had an effect on the real cost 
mie t OT). 

These new parameter values will be fed into’ the 


controller and function minimization subroutine as well. 
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Figure 6.1 Adaptive Control Scheme 


C. ADVANTAGES 


* Using this adaptive scheme we don't use sensors which 


may be unrealistic and the use of predetermined filters 


which are so expensive can be avoided. 
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Figure 6.2 Decision Generator 


¢ Provides a good choice for the function minimization 
Subroutine to start it's iteration algorithm and there- 


fore we get the optimum values as quickly as possible. 


D. DISADVANTAGES 


* The search algorithm must be carefully determined. 
¢ Since we don't include sensor's measurements there is 


no indications on how to perform the search. This 
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problem can be eliminated when the NAVSTAR/GLOBAL 
POSITION SYSTEM (GPS) will be able to provide precision 


navigation data. 
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A. 


VIL. CONCLUSIONS AND RECOMMENDATIONS 


CONCLUSIONS 


A properly designed controller can minimize the rudder 
activity providing desired overseas heading as well,and 
therefore does result in substantial improvement in 
propulsion efficiency. 

Actual savings in fuel cannot yet be determined since 
there is no information available from the conventional 
autopilot and therefore there is no possibility for 
comparison. 

It is believed that the performance index used in this 
research is a fairly adequate function.Doubts arise 
from the weighting factor which is included and this 
because lambda is based on assumptions and it's accu- 
racy 1S not certain. 

BDeeddaptsve controller which minimizes propulsion 
losses due to steering is needed when environmental 
conditions and ship characteristics change. 

Studying all the investigated sailing modes, it turns 
out that the cost surface is flat.As a consequence 
determination of the controller parameters does not 
require high accuracy. 

The study shows that the use of the third order ship 
Nomoto model is a reasonable choice instead of using 
the ship's equation of motion,which involves both the 
Sway and yaw equations. 

It can be assumed from the work done in this research 
that some of the findings could have applicability to a 
number of ship types.However it is not possible to make 


inference, of a general nature concerning all ships 
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De 


from the results of the SL-/ which represents a partic- 


ular type of a particular class Geeshitps. 


RECOMMENDATIONS FOR FUTURE STUDY 


Additional work has to be done for obtaining optimum 


controller parameters under an expanded range of oper- 


ating conditions. The more optimum controller param- 
eter values available, the better the determination of 
the search algorithm recommended in the adaptive 


control scheme. 


e A study including the surge equation in our ship model 


is recommended since in reality added resistance due to 
steering reduces ship's speed and so far we assumed 
constant speed. In addition with the use of the surge 
equation we can determine actual energy losses. 
Investigation of more advanced methods of adaptive 
control based on "on-line" determination of plane 
characteristics. 

So far we were interested in minimizing rudder and 
yawing activity to reduce propulsion losses using a 
particular performance Criterion in open 
seas.Considering some other types of control such as 
automatic control for replenishment requires an inves- 
tigation as to the suitability of the cost function, 


and a comparison with other potential cost functions. 
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APPENDIX A 


DETERMINATION OF NOMOTO THIRD ORDER MODEL-BOXPLX 


//PROGRA JOB (????,0356), 'RESEARCH' ,CLASS=G 


// EXEC FORTXCG, PARM.FORT='OPT(2)',IMSL=DP, REGION=1024K 
//FORT.SYSIN DD * 

C THIS PROGRAM WILL OBTAIN THE CONTROLLER OPTIMAL GAINS. 

C IT IS REFERENCED IN CHAPTER 5. 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
C OBTAINED CHANGE XS(*) TO X(*) AND DELETE XU(*),AND XL(*). 


fo XS (1) 
Gee XL(1) 
eee xU(I) 


DIMENSION XS(5),XU(5),XL(5) 
XS(1)=1.7466650 

XS(2)=35.2978973 

XS(3)=22.2485657 

XS(4)=13.6511078 

XS(5)=22.1172638 

IS THE STARTING GUESS 

IS THE LOWER LIMIT FOR THE I'TH VARIABLE 
IS THE UPPER LIMIT FOR THE I'TH VARIABLE 
XL(1)=0.1 

XU(1)=4.0 

XL(2)=0.1 

XU(2)=80.0 

XL(3)=1.0 

XU(3)=50.0 

XL(4)=1.0 

XU(4)=30.0 

XL(5)=1.0 

XU(5)=50.0 


C A DESCRIPTION OF THE FOLLOWING PARAMETERS 
C IS DISCUSSED IN BOXPLX 
R= OS: 


65 


NTA=1000 
NPR=100 
NAV=0 
NV=5 
IP=0 
C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
C CALL PLANT(X) 
C IF ONLY SIMULATION IS WANTED 
CALL BOXPLX(NV,NAV,NPR,NTA,R,XS,IP,XU, XL, YMN, IER) 
WRITE (6,25) 


25 FORMAT(1X,' OPTIMAL GAINS"',/) 
DO 30 I=1,5 

30 WRITE(6,40)1,XS(L) 

40 FORMAT(1X, K( . 1226 = Sr ae 
STOP 
END 


SUBROUTINE PLANT(XX) 

C SUBROUTINE PLANT(XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL" 8 Lu le2b3 4 ome 
REAL*8 X,XDOT,Y,YDOT,U,UDOT,V, VDOT, YAW,R,RDOT 
REAL*8 TIME, ETIME, XUDOT, XUU,XVR,XVV,XDD 
REAL*8 YV,YR,YD,YVVR,YVRR,YVVV,YRRR, YDDD, YVDOT 
REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
REAL*8 RHO,IZ,FX,FY,MZ,XP,MASS, DELT,MZI,RXI,WA,WE 
REAL*8 DYAWE, YAWE, YAWC,ISE,ISR,LAMDA,D,RYR,RYI,MZR 
REAL*8 K1,T1,T2,D,X2,DX2,S,RX,RY,RZ,TX,TY,TZ,RXR 
REAL*8 T3,T4,X3,DX3,X4 
DIMENSION XX(5) 


CLOSE LOOP ANALYSIS WITH FILTER 
INITIAL CONDITIONS FOR INTEGRATION 


SIMULATION END TIME IN SECONDS 
ETIME=600 .0 


Gii''G@)- Gi ©) C3 
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C 


C 


C 


TIME=0.0 
LeOunm al 
INITIALIZE THE COST FUNCTION 
ISE=0.0 
ISR=0.0 
ED OO 
LAMDA=8.128 
GAIN COEFFICIENTS TO BE OPTIMIZED 
l= xed) 
re Kee) 
T2= KS) 
T3=XX(4) 
a= kN (5) | 
WRITE(6,1010) K1,T1,T2 


C1010 FORMAT(1X,'Kl =',F15.7,'Tl =',F15.7,' T2 =',F15.7) 


C 


X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
0 
XDOT=0.0 
YDOT=0.0 

U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 
RDOT=0.0 

ORDERED SPEED IN FEET/SEC 

38.82 FI/SEC=23 KNOTS 

UG=38.82 

AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
eae 

D = RUDDER ANGLE 
D=0.0 
L=880.5 


6/ 


C 


L2=L**2 
L3=L*L*L 
aie 
L5=L*L4 
L6=L*L5 


SEA DISTURBANCE 


C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
C MOMENTS IN Z . 


Clim Giencae co OC) Cy CC) OC) (Cn) Ci wil OG) VO) Cy Oe) 


FX=0. 

FY=0. 

MZ=0. 
RXR=-0.15744D+05 
RXI=-0.19950D+06 
RYR=0.52365D+04 
RYI=0.18699D+06 
MZR=-0.29870D+08 
MZ1=— 0) 3575 0D+07 
RXR=-0.50230D+04 
RXI=0.12712D+05 
RYR=0.35290D+04 
RYI--0.31909D+05 
MZR=0.38826D+07 
MZI=-0.64313D+07 
RXR=0.28540D+04 
RXI=-0.99574D+04 
RYR=-0.85441D+04 
RYL-0).. 39595D205 
MZR=-0.13014D+08 
MZI=0.11348D+08 
RXR=-0.75642D+04 
RXI=0.83497D+04 
RYR=0.23379D+05 
RYI=-0.81502D+05 
MZR=0.28622D+07 
MZI=-0.19388D+08 
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C 


RXR=-0.37916D+04 
RXI=0.16381D+04 
RYR=-0.76647D+05 
RYI=-0.37685D+05 
MZR=-0.83915D+07 
MZI=-0.53176D+07 
RX=DSQRT (RXR**2+RXI**2 ) 
RY=DSQRT(RYR**2+RY1I**2 ) 
RZ=DSOQRT (MZR*¥*2+MZ1**2 ) 
TX=DATAN(RXI/RXR) 
TY=DATAN(RYI/RYR) 
TZ=DATAN(MZI/MZR) 
SIGNIFICANT WAVE HEIGHT; SEA STATE 1-5,2-10,3-15 
feb 22.5 46-27 ,7-35,8-42,9-60 
WA=17.5 
ENCOUNTER FREQUENCY .1,.2,.3,.4,.5,-6,-75,1.0,1.5,2.5 
WE=0.2 
HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS 
PARAMETERS 
RHO=1.9876 
MASS=(.0044)*(.5*RHO*L3) 
IZ=(0.00028 )*(.5*RHO*LS ) 
YAWE=0.0 
X2=0.0 
DX2=0.0 
5S 100 
pres) (0 


200 CONTINUE 


S=DSQRT (U**2+V¥*2 ) 
INPUT YAW COMMAND 
YAWC=0.0 
IF (TIME.GE.0.0) YAWC=0.0 


C ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL - YAW ORDERED ) 


C 


( COMPENSATOR FILTER ) 
YAWE=YAW - YAWC 
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C 
C 


CP.1' Gy; On Cl, 462! 6) 


qdNaany 


DX2=(YAWE-X2)/T2 
X4=K1*(T1*DK2+X2) 
DX3=(X4-X3)/T4 
D= (T3*DX3+X4) 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE , SPEED , | 
ENCOUNTER FREQUENCY 
XUDOT= (- .0001)*(.5*RHO*L3) 
KU= (-0.0253)*(.5*RHO*1L2*S) 
XUU= (-0.0003)*(.5*RHO*L2 ) 
XVR=(0.0039)*(.5*RHO*L3) 
XVV=(-.0012)*(.5*RHO*L2) 
XDD= (-0.0005)*(.5*RHO*L2*S**2 ) 
LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758)*(.5*RHO*L2*S) 
YR= (0.0023)*(.5*RHO*L3*S) 
YD=(0.00145)*(. 5*RHO*L2*S**2 ) 
YVVR=(0.01)*(.5*RHO*L3/S) 
YVRR=(-0.008)*(.5*RHO*L4/S) 
YVVV=(-0.03)*(.5*RHO*L2/S) 
YRRR=(0.003)*(.5*RHO*L5/S) 
YDDD= (-0.0005)*(.5*RHO*L2*S**2) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE 
CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 
ENCOUNTER FREQUENCY 
YVDOT=(-0.0039)*(.5*RHO*L3) 
SPEED=23 KNOTS, ENCOUNTER ANGLE = 30 , ENCOUNTER 
FREQUENCY =0.2 
YVDOT=-0.30908D+07 
YV=-0.81271D+04 
YVDOT=-0.36185D+07 
YV=-0.24757D+06 
YVDOT=-0.32890D+07 
YV=-0.11775D+07 
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YVDOT=-0.23038D+07 
YV=-0.18267D+07 
YVDOT=-0.59800D+06 
YV=-0.13260D+07 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO*L3*S) 
NR=(-0.00105 )**(.5*RHO*L4*S ) 
ND=(-0.0007)*(.5**RHO*L3*S**2 ) 
NVVR=(-0.015)*(.5**RHO*L4/S) 
NVRR=(-0.008)*(.5*RHO*L5/S) 
NVVV=(0.01)*(.5*RHO*“L3/S) 
NRRR=(-0.006)*(.5%*RHO*L6/S) 
NDDD=(0.0001)*(.5*RHO*L3*S**2 ) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE 
CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 
ENCOUNTER FREQUENCY 


NRDOT=(-0.00027)*(.5*RHO“LS ) 
SPEED=23 KNOTS, ENCOUNTER ANGLE = 90 , ENCOUNTER 
FREQUENCY =0.2 

NRDOT=-0.26251D+12 

NR=-0.53637D+09 

NRDOT=-0.20125D+12 

NR=-0.94970D+10 

NRDOT=-0.18671D+12 

NR=-0.46860D+11 

NRDOT=-0.14518D+12 

NR=-0.87538D+11 

NRDOT=-0.37261D+11 

NR=-0.69856D+11 
REGULAR WAVE SEA STATE 

FX=WA*RX*DCOS (WE*TIME+ TX ) 

FY=WA*RY“DCOS (WE*TIME+TY ) 

MZ=WA*RZ*DCOS (WE*TIME+TZ) 


C U ACTUAL SPEED 


vl 


UC COMMANDED SPEED 


XP 


= PROPELLER THRUST 
XP=- XUU*UC**2 


EQUATIONS OF MOTION 


UDOT=( (KVR + MASS)*V*R + XUU*U**2 + XVV*V2"2 

1 + XDD*D*D + FX + XP )/(MASS-XUDOT) 

VDOT=(YV*V + (YR-MASS*U)*R + YD*D + YVVR*V**2*R 
1 +YVRR*V*R**2 +YRRR*R**3 

2 + YDDD*D**3 + FY )/(MASS-YVDOT) 

RDOT=( NV*V + NR*¥R + ND*D + NVVR*V**2"R 

1 + NVRR*V*R#*2 +NVVV"V%23 

2 + NRRR*R**3 + NDDD*D**3 + MZ )/(IZ-NRDOT) 


WHEN TO PRINTOUT 


IF (ICOUNT. EO: Lip mcomro) 56 
GO TO 300 


CONVERT RADIANS TO DEGREES 


50 


YAWDEG= YAW*57.296 

RDEG=R*57.296 

RDDEG=RDOT*57 .296 

DDEG=D*57.296 

YAWC=YAWC*57 .296 

WRITE (6,100) TIME,XP,X,XDOT,Y,YDOT 

1 ,UC,U,UDOT,V,VDOT, YAWC, YAWDEG , RDEG , RDDEG , DDEG 


100 FORMAT(1X, TIME=' ,F8.375 "SEC ™Xe— — elie2. J oee 


lL X=',F8.2, 

1’ FT xXDOT=',F8.4,' FT/SEC Y=',F8.2 Rte evwDen: 
leans ot, SEG 

1',/,2X,' UC=',F8.4,° FI/SEC U=" heeceee one 
1 UDOT=',F10.6, 

1 * FI/SEC®*2  V=',F8.4,' FI/SEG ME DOT— lamer 
lL 2B SEG 

1 ,/,2X,'YAWC=',F8.4,' DEG YAW=’ 

1  ,F15.7,' DEG YAW RATE=',F15.7,' DEG/SEC 

1 YAW ACCEL=',F15.7,' 

1 DEG/SEC**2',/,2X,'RUDDER =',F15.7,' DEG ',/) 
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ICOUNT=1 
C TEST IF WANT TO STOP 
300 IF (TIME.GE.ETIME) GO TO 400 
C INTEGRATION STEP SIZE DELT 
Dee —eo 
C INTEGRATION 
U=U+UDOT*“DELT 
V=V+VDOT*DELT 
R=R+RDOT*“DELT 
YAW=YAW+R*DELT 
¥2=X2+DX2*DELT 
X3=X3+DX3*DELT 
C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
C XDOT=U*DCOS (YAW) -V*DSIN (YAW) 
C YDOT=U*DSIN(YAW)+V*DCOS (YAW) 
C X=X+XDOT*“DELT 
c Y=Y+YDOT*DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D**2 
GO TO 200 
G J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 
WRITE(6,500) ISE,ISR,TDIFF,K1,T1,1T2,T3,T4 
500 FORMAT(’ ',1X,'TOTAL=',F15.7,2x, 
Halas 7.28. °T1=",F15.7,2X%,'T2=',F15.7,2X; 
leet amewe 5 7 Oke TAS P15 27.) 
RETURN 
END 
SUBROUTINE BOXPLX (NV,NAV,NPR,NTZ,RZ,XS,IP,BU,BL, 
LYMN , IER) 


DIMENSION V(50,50), FUN(50), SUM(25), CEN(25), 
1XS(NV),BU(NV) , BL(NV) 
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KV = 5 


EP = 1.E-6 

NTA = 2000 

IF (NTZ.GT.0) NTA = NTZ 
R = RZ 


IF (R.LE.0. GRO RsCE eee ey oe 
NVT = NV+NAV 


TOTAL VARS, EXPLICIT PLUS IMPLICIT 
NT = 0 
CURRENT TRIAL NO. 


NPT = O 
CURRENT NO. OF PERMISSIBLE TRIALS 
NTFS = 0 


CURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGED 


CHECK FEASIBILITY OF START POINT 


DO 4 I=1,NV 


VT = XS(I) 
TRe(BL( 1) SLE. Vi cOmromsl 
iil 
VT = BL(I) 
GO TO 

1 IF (BU(L)=GE.VD)Go Tome 
II =I 
VT = BUCH) 


2 IF (NPR.GT.0O) WRITE (6,49) II 
30 Gee eal 
CEN(I) = VT 
TE (LE bem mcOnEon4 
BL(I) = BL(1I)+AMAX1(EP,EP*ABS(BL(I) )) 
BU(I) = BU(I)-AMAX1(EP,EP*ABS(BU(I))) 
4 SUM(1) = VT 
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NCE 
NUMBER OF CONSTRAINT EVALUATIONS 
a eL 
TF (KEC ler) ) .nOr0)) GO TO 5 
IF (NPR.LE.0) GO TO 12 
WRITE (6,50) 
GO hom 
5 NFE = 1 


NUMBER OF VERTICES (K) = 2 TIMES NO. OF VARIABLES. 
K = 2*NV 


NUMBER OF DISPLACEMENTS ALLOWED. 
NLIM = 5*NV+10 


NUMBER OF CONSECUTIVE TRIALS WITH UNCHANGED 
FE TO TERMINATE. 

NCT = NLIM+NV 

ALPHA = 1.3 

FK = K 

FKM = FK-1. 

BETA = ALPHA+1. 


INSURE SEED OF RANDOM NUMBER GENERATOR IS ODD. 
IQR = R*1.E7 
IF (MOD(IQR,2).EQ.0) IQR=IQR+101 


SET UP INITIAL VERTICES 
FUNG nye= FE(VCL,1)) 
YMN = FUN(1) 
6 FI = 1. 
FUNOLD = FUN(1) 


i 


QO 


©.) 


/ 


DOM L=29e 


Hie EF Le 
LIMT = 0 
LIMT = LIMEes 


END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 


IF (LIMT. GE) NEG colonel 


DO 8 J=1,NV 


RANDOM NUMBER GENERATOR’ (RANDU) 


LO 


Ll 


EZ, 


IQR = IQR*65539 

IF (IQR.LT.O) IQR = IQR+2147483647+1 
RQX = IQR 

RQX = RQX*.4656613E-9 

V(J,1) = BL(J)+RQX*(BU(J)-BL(J)) 

IF (IP.EQ.1) V(J,I)=AINT(V(J,1I)+.5) 
CONTINUE 


DO 10 L=1,NLIM 
NCE = NCE+1 
TE (REC Le) Ee OuOn mC OmmOmanS 


DO 9 J=1,NV 

VT = .5*(V(J,1)+CEN(J)) 

IF (IP.EQ.1) VT = AINT(VT+.5) 
V(J,1) = VT 

CONTINUE 


CONTINUE 


ENE REE 0 \iCOrmO nado 

WRITE (6,51) I 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,1,FUN,CEN,I) 
IER = -l 
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GO TO 48 


ise DO laed— iy NV 
SUM(J) = SUM(J)+V(J,I1) 
14 CEN(J) = SUM(J)/FI 


TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 
NCE = NCE+1 
IF (KE(CEN).EQ.0) GO TO 60 
SUM(J) = SUM(J) -V(J,I) 
GO ale 7 
60 NFE = NFE+1 
FUN(IL) = FE(V(1,I1)) 
15 CONTINUE 


END OF LOOP SETTING OF INITIAL COMPLEX. 
IF (NPR.LE.O) GO TO 17 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K, FUN,CEN,0) 


PiNb> THE WORST VERTEX, THE J THe 
J; = al 


DO 16 [T=2,K 
IF (FUN(J).GE.FUN(I)) GO TO 16 
oi 

16 CONTINUE 


BASIC LOOP. ELIMINATE EACH WORST VERTEX 
IN TURN. IT MUST BECOME NO LONGER WORST,NOT 
MERELY IMPROVED. FIND NEXT-TO-WORST VERTEX, 
THE 'JN'TH ONE. 
17 JN = 1 

IF (J.EQ.1) JN = 2 


DO 18 I-1,K 
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[F (1. EO.g cen ToOmls 
IF (FUN(JN).GE.FUN(I)) GO TO 18 
JIN = 1 

18 CONTINUE 


LIMT = NUMBER OF MOVES DURING THIS TRIAL TOWARD 
THE CENTRIOD DUES TOC SFUNGTEONS YALE 
eM eal 


COMPUTE CENTROID AND OVER REFLECT WORST VERTEX. 


DO 19 I=1,NV 

Viv Vel 

SUM(I) = SUM(1)-VT 

CEN(I) = SUM(I)/FKM 

VT = BETA*CEN(I)-ALPHA*VT 

IF (IP.EQ.1) VT = AINT(VT+.5) 


INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 
19 V(I,J) = AMAX1(AMIN1(VT,BU(I)),BL(I)) 


Nis= Niet 
CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 


20) DORZ aN NEE 
NCE = NCE+1 
EF (KE(CV (1) -EOe0 ce ie e2c 


EVERY 'KV'TH TIME, OVER-REFLECT THE OFFENDING 
VERTEX THROUGH THE BEST VERTEX. 

IF (MOD(N,KV).NE.0) GO TO 22 

CALL FBV (K,FUN,M) 


DO 21 I=1,NV 
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VT = BETA*“V(1,M)-ALPHA*V(I,J) 
IF (IP.EQ.1) VT = AINT(VT+.5) 


21 V(I,J) = AMAX1(AMIN1(VT,BU(I)),BL(I)) 


GO TO 24 


CONSTRAINT VIOLATION: MOVE NEW POINT TOWARD CENTROID. 


22 DO 23 I=1,NV 


VT = .5*(CEN(L)+V(I,J)) 
IF (IP.EQ.1) VT = AINT(VT+.5) 
7 Clee lee idl 


23 CONTINUE 


24 NT = NT+1l 
25 CONTINUE 


heen 


CANNOT GET FEASIBLE VERTEX BY MOVING TOWARD CENTROID, 
OR BY OVER-REFLECTING THRU THE BEST VERTEX. 
EERE O)) cor TO 42 
VRUDEG6m>2 )) NT. J 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,J) 
GOnTOn 42 


FEASIBLE VERTEX FOUND,EVALUATE THE OBJECTIVE FUNCTION. 


26 NFE = NFE+1 


FUNTRY = FE(V(1,J)) 


TEST TO SEE IF FUNCTION VALUE HAS NOT CHANGED. 
AFO = ABS(FUNTRY-FUNOLD) 
AMX = AMAX1(ABS(EP*FUNOLD) , EP) 


ACTIVATE THE FOLLOWING TWO STATEMENTS 
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FOR DIAGNOSTICS PURPOSES ONLY. 


99 


aa 


WRITE (6,99) J,AFO,AMX,FUNTRY,FUNOLD,FUN(J), 


LFUN(JN) ,NTFS,N 


FORMAT (1X,13,6E15.7,215) 

IF (AFO.GT.AMX) GO TO 27 

NTFS = NTFS+1 

IF (NEFS2LT NGI) meommon 2c 

IER = 0 

IF (NPR.LE.0O) GO TO 42 

WRITE (6,53) K 

CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,0) 
G0 Tom 

NTFS = 0 


IS THE NEW VERTEX NO LONGER WORST? 


28 


IF (FUNTRY.LT.FUN(JN)) GO TO 34 


TRIAL VERTEX IS STILL WORST; ADJUST TOWARD CENTROID. 
EVERY ‘KV'TH TIME, OVER-REFLECT THE OFFENDING 
VERTEX THROUGH THE BEST VERTEX. 


Za 


30 


LIMT = LIMT+1 
IF (MOD(LIMT,KV).NE.0) GO TO 30 
CALL FBV (K,FUN,M) 


DO 29 I=1,NV 

VT = BETA*V(1I,M)-ALPHA*V(I,J) 

iF (LP 20s), Vial inG@a oe 

V(I,J) = AMAX1(AMIN1(VT,BU(IL)),BL(I1)) 


GO TO 32 


DO 31 I=1,NV 

VT = .5*(CEN(1)+V(I,J)) 

TF (IP £Onl) VI = AINGC Tae 
AGL) = wit 
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31 CONTINUE 


2 1k (LoMm. kT. NLEM) Go. tO. 33 


CANNOT MAKE THE 'J'TH VERTEX NO LONGER WORST 

BY DISPLACING TOWARD OVER-REFLECTING 

THRU THE BEST VERTEX. 
epee 
AGN Ree EPO CON TOn42 
WRITE (6,52) NT, J 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN, J) 
CO WO AP 

33 NT = NT+1 
Counon2@ 


SUCCESS: WE HAVE A REPLACEMENT FOR VERTEX J. 
34 FUN(J) = FUNTRY 

FUNOLD = FUNTRY 

NPT = NPT+1 


Even LOO TH PERMISSIBLE TRIAL, RECOMPUTE 
CENTRIOD SUMMATION TO AVOID CREEPING ERROR. 


IF (MOD(NPT,100).NE.0) GO TO 37 


DO 36 I=1,NV 


SUM(I) = 0. 
Bommel k 
35 SUM(I) = SUM(I)+V(I,N) 


CEN(L) 
36 CONTINUE 


SUM(1)/FK 


LC = Q 
COPTOrs 7 
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OQ 


37 DO 38 I=1,NV 
38 SUM(I) = SUM(1)+V(I,J) 


LC = J 


39 IF (NPR.LE.0) GO TO 40 
IF (MOD(NPT,NPR).NE.0) GO TO 40 


CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,CEN,LC) 
HAS THE MAX. NUMBER OF TRIALS BEEN REACHED 
WITHOUT CONVERGENCE? 

IF NOT, GO TO NEW TRIAL. 


40 IF (NT.GE.NTA) GO TO 41 


NEXT-TO-WORST VERTEX NOW BECOMES WORST. 


J = JN 
GO TO 17 
41 IER = 3 


IF (NPR.GT.0) WRITE (6,54) 


COLLECTOR POINT FOR ALL ENDINGS. 
1) CANNOT DEVELOP FEASIBLE VERTEX. IER = 


iE 
2) CANNOT DEVELOP A NO-LONGER-WORST VERTEX.IER = 2 
3) FUNCTION VALUE UNCHANGED FOR K TRIALS. IER = 0 
4) LIMIT ON TRIALS REACHED. IER = 3 
5) CANNOT FIND FEASIBLE VERTEX AT START. TER joes I 


42 CONTINUE 


FIND BEST VERTEX. 
CALL FBV (K,FUN,M) 
IF (IER.GE.3) GO TO 44 


RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER 
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THAN THE PREVIOUS, OR IF THIS IS THE FIRST TRY. 
IF (NPR.LE.O) GO TO 43 
WRITE (6,55) (M,YMN,FUN(M)) 
43 IF (FUN(M).GE.YMN) GO TO 47 
IF (ABS(FUN(M)-YMN).LE.AMAX1(EP,EP*YMN)) GO TO 47 


GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED. 
44 YMN = FUN(M) 
FUN(1) = FUN(M) 


DO 45 I=1,NV 

CEN(I) = V(I,M) 

SUM(I) = V(I,M) 
45 V(I,1) = V(I,M) 


DO 46 I=1,NVT 
46 XS(I) = V(I,M) 


LUeGlERGE 3) CO TO 6 

47 IF (NPR.LE.O) GO TO 48 
CALL BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FUN,V(1,M),-1) 
WRITE (6,56) FUN(M) 

48 RETURN 


49 FORMAT (SOHOINDEX AND DIRECTION OF 
LOUTLYING VARIABLE AT STARTIS) 

50 FORMAT (S5OHOIMPLICIT CONSTRAINT 
LVIOLATED AT START. DEAD END. ) 

51 FORMAT ('OCANNOT FIND FEASIBLE’ ,I4,'TH 
1IVERTEX OR CENTROID AT START.) 

52 FORMAT (1LOHOAT TRIAL 14,54H CANNOT FIND 
LFEASIBLE VERTEX WHICH IS NO LONGER 
1WORST,1I4,15X, 'RESTART FROM BEST VERTEX. ') 

53 FORMAT (40HOFUNCTION HAS BEEN ALMOST 
LUNCHANGED FOR 15,7H TRIALS) 
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54 FORMAT (27HOLIMIT ON TRIALS EXCEEDED. ) 
55 FORMAT ('OBEST VERTEX IS NO.',13,' 
1 OLD MIN WAS 'E15.7', NEW MIN IS ',E15.7) 
56 FORMAT ('OMIN OBJECTIVE FUNCTION IS ',E15.7) 
END 
SUBROUTINE FBV (K,FUN,M) 
DIMENSION FUN(50) 
M = 1 


DO 1 I=2,K 
IF (FUN(M).LE.FUN(I)) GO TO 1 
M=I 

1 CONTINUE 


RETURN 

END 

SUBROUTINE BOUT (NT,NPT,NFE,NCE,NV,NVT,V,K,FN,C,IK) 
DIMENSION V(50,50), FN(50), C(25) 

WRITE (6,4) NT,NPT,NFE,NCE 


DO 1 I=1eK 
WRITE (6.5). EN(L). (Viable 
IF (NVT.LE.NV) GO TO 1 
NVP = NV+l 
WRITE (6,6) (V(J,1),J=NVP,NVT) 
1 CONTINUE 


iF (1k. NERO) eco noe? 


WRITE (6,79 (CULL) 1-12 in 
RETURN 

9 TE ((DKACE 0) aCe on 3 
WRITE (6,8) (C(I),I=1,NV) 
RETURN 

3 WRITE (6,9) IK,(C(1I),I=1,NV) 
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RETURN 


eEORMARMGZONO Onan mincm = 915, 4X. 
L'NO. FEASIBLE TRIALS = ',15,4X,'NO. FUNCTION 
LEVALUATIONS = ',15,4X,'NO. CONSTRAINT EVALUATI 
WONS = 15/0 FUNCTION VALUE' ,6X, ' INDEPENDENT 
1VARIABLES/DEPENDENT OR IMPLICIT CONSTRAINTS' ) 
RORMAT (IHR IS) 7eoN 7 bl’. 7/(21X,7EI4.7)) 
FORMAT (21X,7E14.7) 
FORMAT (1OHOCENTROID 11X,7E14.7/(21X,7E14.7)) 
FORMAT ('O BEST VERTEX' ,7X,7E14.7/(21X,7E14.7)) 
FORMAT ('OCENTROID LESS VX',12,2X,7E14.7/ 
(DIG nO 
END 

FUNCTION FE(X) 

DIMENSION X(5) 

COMMON TDIFF 

CALL PLANT(X) 

FE=TDIFF 

RETURN 

END 

FUNCTION KE(X) 

DIMENSION X(5) 

KE=0 

RETURN 

END 
//GO.SYSIN DD * 
/* 


Oo ON HD 
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A. 


CONTROELERS JG. 





Figure B.1l 





APPENDIX B 
CONTROLLER DESIGN 


Ty 
i + 1 
a eee +. 
T3 -£ S 
i 
T3 


Block Diagram of Controller C 
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Pepuice Beleeecorrespotids to controller ‘CC’ which has the 


Grin : 


K,(1+T 7S) (1+TyS) 





CT) Clie Ss) 
CBr) 
Verifying that equation B.1l corresponds to controller 


me swe have: 








K 
Bs (B.2) 
E +K74T 7D 
x - 
a - DX9=SX9 (B.3) 
x 
oe : AQ = 
X¥, itt2s 1+7T25 
(B.4) 


Substituting equations B.3 , B.4 into equation B.2 we 


have: 





XK Kt a 
A= ——_—_ +  ____..__ 
iS i>. 
(B53) 
We also have: 
A : 
= + T,, DX 
eee S ime (B.6) 
X 1 = 
ao Go DX 3=%39 (B.7) 
DX 3 5 
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i (B.8) 





Substituting equations B.7 , B.8 into equation B.6 we 
have: 


A, Teas 


SU Teese SURE 
(B. 9} 


Now substituting equation B.5 into equation B.9 we 


have: 
Poy A A ies i 
~(1+T5S)(1+#T3S) (1+#T2S)(1+T3S) 
X1K,T,$ X1K1T1TyS 
+ 


Finally rearranging terms equation B.10 can be 


written: 


Xy Ky (1+TySt+T,S+TyT 8°) Ky (14TyS) (147, S) 
A Gees One) (lst 5S) Claes (B21ap) 


B. CONTROLLER "B"” 


Figure B.2 corresponds to controller 'B' which has the 


form: 
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Figure B.2 Block Diagram of Controller B 


K,(1+T 5S) 
helio o) (lt T3s) 


(B.12) 


Verifying that equation B.12 corresponds to Controller 


"B' we have: 


KX 
Az-l°1l , 
a (B.13) 
be 
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SS DX9= X98 (B.14) 
DXy 9° 
Xo x1 
1 
——= X 2 2-—aee (B. lay 
X, T2Stl 1+#7T)8 
Substituting eguations B.14 , B.15 into equations 
we have: 
- a ae 
1+T,S 1+T,S (B.16) 


In a similar way we also can derive equation B.1/ 





14+T38 
(B.17) 


= 


Substituting equation B.16 into equation B.1/ we have: 


XQ sate ee 
" Ceme Geiss) | (Genser. (B.18) 


Finally rearranging terms equation B.18 becomes: 


A eee 


¥ (1#T78)(1+T38) (Pa 
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fee CONTROLLER “A” 





Figure B.3 Block Diagram of Controller A 


Figure B.3 corresponds to controller ‘A’ which has the 


form: 


i Ss) 


1+T,S 
(B20) 
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Verifying that equation B.20 corresponds to Controller 


'A' we have: 





= ae + KT 1DX 
AG 1+T)5 141449 
(B. 2a) 
But since we know that: 
__ a ee DX9=X9S 
x2 5S (B.22) 
x 
ome a 


SubStltuting equatiom Bacce. B.23°into equation Baas 





we have: 
ae XK, . K,T,%45 
as ies (B.24) 


Finally rearranging terms equation B.24 becomes: 


pc aariee 


X, its (B.25) 


ele 


CODING OF THE EQUATIONS 


= 2 = Por Controllers © -s22—--------- 





YAWE = YAWC - YAW 
YAWE = YAW - YAWC 
(YAWE - Xo) 
Dxq = 
i 


Xo oa Xo + DX>° DE eT 
j= Ky CX + T,DX5) 


A-X3 X3 
T3 


X3 + DX 3° DE 


il 


io 

~ 
Ow 
i 


i = B3 + T, Dx 


2255 For Controller "B"------------- 


Ipeecracion-------= 
YAWE = ea ee YAWC 
YAWE ~ x ; 
XQ = ———____° = Xo + DX> = DEET 
T2 
Neg a ot Do) = ies 
A= Ki(X2 + TDX>) 
D = 
xX 
— Bomm@omeroller “A"--2---4--4---Integration-------< 
YAWE = YAW - YAWC 
YAWE - x, 
Im eee eD Kerr DEL 
T9 


Do s= Kz (X94 + i DX) 


93 


For all of the above cases the equations include the 
error detector and the controller which are indicated in 


Figure B.4 


CORO hehe ae. 
tpt 
Wen 


Fabel a 





Figure B.4 General Scheme of Control 
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APPENDIX C 
RESPONSE OF THE SYSTEM FOR REGULAR SEAS 


//PROGRA JOB (????,0356),'RESEARCH' ,CLASS=A 


// EXEC FORTXCG, PARM. FORT='OPT(2)',IMSL=DP , REGION=1024K 
//FORT.SYSIN DD * 
C 


HAVE BEEN OBTAINED CHANGE XS(*) TO X(*) AND DELETE 
XU(*),AND XL(*). 
COMMON J 
DIMENSION X(5) 
X(1)=1.8287125 
X(2)=1.1652012 
X(3)=10.5659571 
¥(4)=11.7124157 
X(5)=20.5683947 
@e GALL PLANT (X) 
C IF ONLY SIMULATION IS WANTED 
CALL PLANT(X) 
WRITE (6,25) 
25 FORMAT(1X,' OPTIMAL GAINS',/) 
DO 30 I=1,5 
30 WRITE(6,40)1,X(I) 
40 HOnMom@loan' x (' ym2," )s=* -F14.7) 
WRITE(6,50) J 
50 FORMAT(1X,'J = ',E15.10) 
STOP 
END 
SUBROUTINE PLANT(XX) 
C SUBROUTINE PLANT(XX) SIMULATES THE SHIP 


C 
C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS 
C 
c 
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etc - @ Gi. 


C 


COMMON TDIFF 

REAL“8 Lj l2 soeea eos 

REAL*8 X,XDOT,Y,YDOT,U,UDOT,V,VDOT,YAW,R,RDOT 
REAL*8 TIME, ETIME, XUDOT, XUU, XVR,XVV, XDD 

REAL*“8 YV,YR,YD,YVVR,YVRR,YVVV,YRRR,YDDD, YVDOT 
REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
REAL*8 RHO,IZ,FX,FY,MZ,XP,MASS,DELT,MZI,RX1,WA,WE 
REAL*8 DYAWE,YAWE,YAWC,ISE,ISR,LAMDA,D,RYR,RYI,MZR 
REAL*8 K1,T1,T2,D,X2,DX2,S,RX,RY,RZ,TX,TY,TZ,RXR 
REAL*8 T3,T4,X3,DX3,X4 

DIMENSION XX(5) 


CLOSE LOOP ANALYSIS WITH FILTER 


INITIAL CONDITIONS FOR INTEGRATION 
SIMULATION END TIME IN SECONDS 
ETIME=600. 
TIME=0.0 
ICOUNT=1.0 
INITIALIZE THE COST FUNCTION 
ISE=0.0 
ISR=0.0 
TDIFF=0.0 
LAMDA=8 .128 
GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX(1) 
T1=XX(2) 
T2=XX(3) 
T3=XX (4 ) 
T4=XX(5) 
WRITE(6,1010) K1,T1,T2 


C1010 FORMAT(1X,‘K1 =',F15.7,'Tl =".Fl5ey) 12 =o 7) 
C X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 


X=0.0 
Y=0.0 
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C 


cae) 


Cy <@) 


OQ OQ 467 @e2 ae 


XDOT=0.0 
YDOT=0.0 
U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 
RDOT=0.0 
ORDERED SPEED IN FEET/SEC 
38.82 FT/SEC=23 KNOTS 
UC=38.82 
AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 
D = RUDDER ANGLE 
D=0.0 
L=880.5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*L5 
SEA DISTURBANCE 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
MOMENTS IN Z 
FX=0. 
FY=0. 
MZ=0. 
RXR=-0.15744D+05 
RXI=-0.19950D+06 
RYR=0.52365D+04 
RYI=0.18699D+06 
MZR=-0.29870D+08 
MZI=-0.35751D+07 
RXR=-0.50230D+04 
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© CO) © © '@ 


Gy “A> GO) a-G). "Ol eG OG, OS Gs OO 


RX -00 127 U2 
RYR=0.35290D+04 
RYT=-0-3 197090205 
MZR=0.38826D+07 
MZI=-0.64313D+07 
RXR=0.28540D+04 
RXI=-0.99574D+04 
RYR=-0.85441D+04 
RY 0=0R 3759 5Da0 
MZR=-0.13014D+08 
MZI=0.11348D+08 
RXR=-0.75642D+04 
RXI=0.83497D+04 
RYR=0.23379D+05 
RYI=-0.81502D+05 
MZR=0.28622D+0/ 
MZI=-0.19388D+08 
RXR=-0.37916D+04 
RXI=0.16381D+04 
RYR=-0.76647D+05 
RYI=-0.37685D+05 
MZR=-0.83915D+07 
MZt == Oreo 176D+07 


RX=DSQRT(RXR**2+RKI**2 ) 
RY=DSQRT (RYR**2+RY1I**2 ) 
RZ=DSQRT (MZR**2+MZ1**2 ) 


TX=DATAN(RXI/RXR) 

TY=DATAN(RYI/RYR) 

TZ=DATAN(MZI/MZR) 
SIGNIFICANT WAVE HEIGHT; 


SEA STATE Jame ver 22. 


“= Le 22 6-277 - 35 6 4 


WA=17.5 


ENCOUNTER FREQUENCY .1,.2,.3,.4,-5,.6,-75,1.0,1.5,2.5 


WE=0.6 


HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE 
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Cc 


OQ? 


GG) GE) 


AS PARAMETERS 


200 


RHO=1.9876 
MASS=(.0044)*(.5*RHO*L3 ) 
IZ=(0.00028)*(.5*RHO*LS ) 
YAWE=0.0 

eG 

DX2=0.0 

X3=0.0 

DX3=0.0 

X4=0.0 

CONTINUE 
S=DSQRT(U**2+V**2 ) 


INPUT YAW COMMAND 


YAWC=0.0 
IF (TIME.GE.0.0) YAWC=0.0 


ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL-YAW ORDERED) 


( COMPENSATOR FILTER ) 


YAWE=YAW - YAWC 
DX2= (YAWE-X2)/T2 
X4=K1*(T1*DX2+X2) 
DX3=(X4-X3)/T4 

D= (T3*DX3+X4) 


AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
FREQUENCY 


XUDOT=(- .0001)*(.5*RHO*L3) 

XU= (-0.0253)*(.5*RHO*L2*S ) 
XUU= (-0.0003)*(.5*RHO*L2) 
XVR=(0.0039)*(.5*RHO*L3 ) 
XVV=(- .0012)*(.5*RHO*L2 ) 

XDD= (-0.0005)*(.5*RHO*L2*S**2 ) 


C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 


C 


YV=(-0.00758)*(.5*RHO*L2*S) 
YR=(0.0023)*(.5*RHO*L3*S ) 


22 


Owe Q GO) G GG) C)1-Q 7 


©). GY GQaie 


YD=(0.00145 )*(.5*RHO*L2*S**2 ) 
YVVR= (0.01)*(.5*RHO*L3/S) 
YVRR= (-0.008 )*(.5*RHO*L4/S) 
YVVV=(-0.03)*(.5*RHO*L2/$) 
YRRR= (0.003 )*(.5*RHO*L5/S) 
YDDD=(-0.0005)*(.5*RHO*L2*S**2 ) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE 
CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED 
ENCOUNTER FREQUENCY 
YVDOT=(-0.0039)*(.5*RHO*L3 ) 
SPEED=23 KNOTS, ENCOUNTER ANGLE = 30 , ENCOUNTER 
FREQUENCY =0.4 
YVDOT=-0.30908D+07 
YV=-0.81271D+04 
YVDOT=-0.36185D+07 
YV=-0.24757D+06 
YVDOT=-0.32890D+07 
YV=-0.11775D+07 
YVDOT=-0.23038D+07 
YV=-0.18267D+07 
YVDOT=-0.59800D+06 
YV=-0.13260D+07 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO“L3*S ) 
NR=(-0.00105)*(.5*RHO*L4*S ) 
ND=(-0.0007)*(.5*RHO*L3*S**2 ) 
NVVR= (-0.015)*(.5*RHO*L4/S) 
NVRR= (-0.008)**(.5*RHO*L5/S) 
NVVV=(0.01)*(.5*RHO*“L3/S) 
NRRR= (-0.006)*(.5*RHO*L6/S) 
NDDD= (0.0001)**(. 5*RHO*“L3*S**2 ) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
FREQUENCY 
NRDOT=(-0.00027)*(.5*RHO“LS ) 


3 


100 


Gi G O10 “ee 


Q 


SPEED=23 KNOTS, ENCOUNTER ANGLE = 30 , ENCOUNTER 
FREQUENCY =0.4 
NRDOT=-0.26251D+12 
NR=-0.53637D+09 
NRDOT=-0.20125D+12 
NR=-0.94970D+10 
NRDOT=-0.18671D+12 
NR=-0.46860D+11 
NRDOT=-0.14518D+12 
NR=-0.87538D+11 
NRDOT=-0.37261D+11 
NR=-0.69856D+11 
REGULAR WAVE SEA STATE 
FX=WA*RX*DCOS (WE*TIME+TX) 
FY=WA*RY*DCOS (WE*TIME+TY ) 
MZ=WA*RZ*DCOS (WE*TIME+TZ ) 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=- XUU*“UC**2 
EQUATIONS OF MOTION 
UDOT=( (XVR + MASS)*V*R + XUU*U**2 + XVVeV%2 
1 + XDD*D*D + FX + XP )/(MASS-XUDOT) 
VDOT=(YV*V + (YR-MASS*U)*R + YD*D 
L + YVVR*V**2*R + YVRR*V*R**2 
L + YVVV*V**3 + YRRR*R**3 + YDDD*D**3 
1 + FY )/(MASS-YVDOT) 
RDOT=( NV*V + NR*¥R + ND*®D + NVVR*V**2*R 
Ll + NVRRYV*R**2 
lL + NVVV*V**3 + NRRR*R**3 + NDDD*D**3 
1 + MZ )/(iz-NRDOT) 


G WHEN TO PRINTOUT 


eC UCOUNE. EO -.82 ) GO TO 150 
GO TO 300 


C CONVERT RADIANS TO DEGREES 


Ou 


50 YAWDEG= YAW*57.296 
RDEG=R*57.296 
RDDEG=RDOT*57 . 296 
DDEG=D*57.296 
YAWC=YAWC* 57.296 
WRITE (6,100) TIME, YAWDEG 
G 1 ,UC,U,UDOT,V, VDOT, YAWC, YAWDEG , RDEG,RDDEG, DDEG 
100 FORMAT(1X,F12.8,1X,F12.8) 
1' FT XDOT=",F8.4,° FDR/SEG  Y-" sac nm 
1 YDOT=',F8.4,' FT/SEC 
1',/,2X,' UC=",F8.4,° "FT/SEC deters se 
leh) SEC UDOT=',F10.6, 
' FT/SEC**2  V=',F8.4,' FT/SEC 
VDOT=" RlOmos FL) SEG. 28 
,/,2X,'YAWC=',F8.4,' DEG YAW=' 
SFIS(7," DEG YAW RATE] sense 
DEG/SEC YAW ACCEL= ' 
nS 7 DEC GEG 0.) kw 
RUDDER =' | FI5.7,° DEG 7) 
WRITE (6,101) TIME, DDEG 
C101 FORMAT(1X,F12.8,1X,F12.8) 
ICOUNT=1 
C TEST IF WANT TO STOP 
300 IF (TIME.GE.ETIME) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT=1.0 
C INTEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
YAW=YAW+R*DELT 
X2=X2+DX2*DELT 
X3=X3+DX3*DELT 
C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS (YAW) -V**DSIN (YAW) 


mG) eG GC) aa Cae GairrC)- (Cy CO) Cy); 
ee oe oe 
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YDOT=U*DSIN(YAW)+V*DCOS (YAW) 
X=X+XDOT*DELT 
Y=Y+YDOT*DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D**2 
GO TO 200 
© J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 
(WRITE Cce SOO KShenGhRarDrPr KieT)] 19.13 .T4 
500 FORMAT(' ',1X,'ISE=',F15.7,' ISR=',F15.7,' 
1) LOPAI—e Els7 2x. 
IeKl= silibeey 2X T=" P15. 752k, 'T2='.F1527,2X, 
Ie uke 5.7 ox 14], F15.7) 
RETURN 
END 
//GO.SYSIN DD * 
/* 


103 


APPENDIX D 
DETERMINATION OF OPTIMAL CONTROLLER PARAMETERS FOR IRREGULAR 
SEAS 


//TRIAL1 JOB (1707,0356),'RESEARCH' ,CLASS=C 
//*MAIN ORG=NPGVM1.1707P 
// EXEC FORTXCG, PARM. FORT='OPT(2)', IMSL=DP, REGION=1024K 
//FORT.SYSIN DD * 
C THIS PROGRAM WILL OBTAIN THE CONTROLLER OPTIMAL 
GAINS. IT IS REFERENCED IN CHAPTER 5. 
IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS 
HAVE BEEN OBTAINED CHANGE XS(*) TO X(*) AND 
DELETE XU(*),AND XL(*). 
DIMENSION XS(5),XU(5),XL(5) 
XS(1)=0.655751 
XS(2)=80.5483 
XS(3)=10.74847 
XS(4)=12.9 
XS(5)=45.09 
C XS(I) IS THE STARTING GUESS 
C XL(I) IS THE LOWER LIMIT FOR THE I'TH VARIABLE 
C xU(I) IS THE UPPER LIMIT FOR THE I'TH VARIABLE 
XL(1)=0.1 
XU(1)=2.5 
XL(2)=40.0 
XU(2)=100.0 
XL(3)=0.1 
XU(3)=20.0 
XL(4)=5.0 
XU(4)=80.0 
XL(5)=60.0 
XU(5)=150.0 
C A DESCRIPTION OF THE FOLLOWING PARAMETERS 


(eS ae ee S| 
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C IS DISCUSSED IN BOXPLX 
R=937 13. 
NTA=1000 
NPR=100 


NAV=0 
NV=5 
IP=0 


C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
©) CALL PLANT(X) 
C IF ONLY SIMULATION IS WANTED 
CALL BOXPLX(NV,NAV,NPR,NTA,R,XS,IP,XU,XL, YMN, IER) 


WRITE 


(6,25) 


25 FORMAT(1X,‘ OPTIMAL GAINS',/) 


DO 30 


iado3 


30 WRITE(6,40)1I,XS(I) 
40 HORMAT CIM Nl ia" E14. 7) 


SLO 
END 


SUBROUTINE PLANT(XX) 
C SUBROUTINE PLANT(XX) SIMULATES THE SHIP 


COMMON 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 
REAL*8 


TDIFF 

L,L2,L3,L4,L5,L6 
X,XDOT,Y,YDOT,U,UDOT,V, VDOT, YAW,R, RDOT 
TIME , ETIME, XUDOT , XUU, XVR, XVV, XDD 

YV,YR,YD, YVVR, YVRR,YVVV, YRRR, YDDD, YVDOT 
NV ,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
RHO ,1IZ,FX,FY,MZ,XP,MASS, DELT 

DYAWE , YAWE, YAWC, ISE,ISR,LAMDA,D 

mien? 13.54. D. x2, DX2, x3, DX3,X4,CH(11),S 


DIMENSION XX(5) 


©) GY Geare 


CLOSE LOOP ANALYSIS WITH FILTER 


INITIAL CONDITIONS FOR INTEGRATION 
SIMULATION END TIME IN SECONDS 
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C 


C 


C 


ETIME=600. 
TIME=0.0 
ICOUNT=1 
INITIALIZE THE COST FUNCTION 
fon= One 
ISR=0.0 
TDIFF=0.0 
LAMDA=8.128 
GAIN COEFFICIENTS TO BE OPTIMIZED 
Kl=xx (1) 
len) 
ED = XR) 
T3=XX(4) 
T4=XX(5) 
WREEE(G, LONG merle a 


C1010 FORMAT(1X, Kl = $8152 73 301 SES 


C 
C 


1 <R2e "FF is27 

X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0.0 © 
XDOT=0.0 
YDOT=0.0 

U,UDOT,V,VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 
RDOT=0.0 

ORDERED SSP EER DaUN = PEE, She 

38.82 FT/SEC=23 KNOTS 

UC=38.82 

AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC 

D = RUDDER ANGLE 
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qa 


©) 


D=0.0 
L=880.5 
162 = ioe 
les = 1)" 15 le 
L4=L*L3 
L5=L*L4 
L6=L*L5 


SEA DISTURBANCE 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
MOMENTS IN Z 


FX=0. 
FY=0. 
MZ=0. 


ISEA IS A SWITCH;ISEA=0 (CALM WATER) ISEA=1 (SEA STATE) 


ISEA=1 


HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS 
PARAMETERS 


200 


RHO=1.9876 
MASS=(.0044)*(.5**RHO*L3 ) 
IZ= (0.00028 )*(.5*RHO*“LS ) 
YAWE=0.0 

X2=0.0 

DX2=0.0 

X3=0.0 

DX3=0.0 

X4=0.0 

CONTINUE 

S=DSQRT (U**2+V"*2 ) 


C INPUT YAW COMMAND 


YAWC=0.0 
IF (TIME.GE.0.0) YAWC=0.0 


C ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL- YAW ORDERED) 
( CONTROLLER FILTER ) 


C 


YAWE=YAW - YAWC 
DX2=(YAWE-X2)/T2 
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Qine iar (= 


C 


GO) OQ: -O). Gl G1 


X4=K1*(T1*DX2+X2) 
DX3=(X4-X3)/T4 
D= (T3*DX3+X4) 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 
XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
FREQUENCY | 
XUDOT= (- .0001)*(.5*RHO*L3) 
XU=(-0.0253 )*(. 5*RHO*E2 en) 
XUU= (-0.0003)*(.5*RHO*L2 ) 
XVR= (0.0039 )*(.5*RHO*L3 ) 
XVV=(-.0012)*(.5*RHO*L2 ) 
XDD=(-0.0005)*(.5*RHO*L2*S**2 ) 
LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758 )*(.5*RHO*L2*S ) 
YR= (0.0023)*(.5*RHO*L3%*S) 
YD=(0.00145)*(.5*RHO*L2*S**2 ) 
YVVR=(0.01)*(.5*RHO*L3/S) 
YVRR=(-0.008)*(.5*RHO*L4/S) 
YVVV= (-0.03)*(.5*RHO“L2/S) 
YRRR=(0.003)*(.5*RHO*LS/S) 
YDDD= (-0.0005)*(.5*RHO*L2*S**2 ) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER 
FREQUENCY 
YVDOT=(-0.0039)%*(.5*RHO*L3) 
SPEED=23 KNOTS, ENCOUNTER ANGLE = ; 
ENCOUNTER FREQUENCY =.75 
YVDOT= -2304300.0 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 
NV=(-0.00213)*(.5*RHO*L3*S) 
NR=(-0.00105)*(.5*RHO*L4*S) 
ND=(-0.0007)*(.5*RHO*L3*S**2) 
NVVR=(-0.015)*(.5*RHO*L4/S) 
NVRR=(-0.008)*(.5*RHO*L5/S) 


108 


NVVV=(0.01)*(.5*RHO*L3/S) 
NRRR=(-0.006)*(.5*RHO*L6/S) 
NDDD= (0.0001)*(.5*RHO*L3*S**2 ) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE 
CHANGED FOR DIFFERENT ENCOUNTER ANGLE , SPEED , 
ENCOUNTER FREQUENCY 
NRDOT=(-0.00027)*(.5*RHO*LS) 
SPEED=23 KNOTS, ENCOUNTER ANGLE = 
ENCOUNTER FREQUENCY =.75 
NRDOT=-1.4518E+11 
C SETS SEA STATE TO ZERO 
IF (ISEA.EQ.1) GO TO 30 
FX=0. 
FY=0. 
MZ=0. 
GO TO 35 
C UNIT 12 HAS THE SEA STATE DATA NAMED CH 
C IT MUST BE SYNCHRONIZED BY TIME 
30 READ (12) CH 
FX=CH(3) 
FY=CH(4) 
MZ=CH(8) 
35 CONTINUE 
U ACTUAL SPEED 
C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
XP=- KUU*UC**2 
C EQUATIONS OF MOTION 
UDOT=( (XVR + MASS)*V*R + XUU*U**2 + XVV*¥V282 
C 1+ XDD*D*D + FX + XP )/(MASS-XUDOT) 
VDOT=(Y¥V*V + (YR-MASS*U)*R + YD*D + YVVR*V*2*R 
1 + YVRR*V*RE*2 
1 + YVVVEV**3 + YRRR*R**3 + YDDD*D**3 
1 + FY )/(MASS-YVDOT) 
RDOT=( NV*V + NR*R + ND*D + NVVR*®V**2*R 


CO)" 26) (Ga CO 


Q 
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L + NVRR*V*R**2 
L + NVVV*V**3 + NRRR*R**3 + NDDD*D**3 
1 + MZ )/(IZ-NRDOT) 

C WHEN TO PRINTOUT 
IF (ICOUNT.EQ.11) GO TO 50 
GO TO 300 

C CONVERT RADIANS TO DEGREES 

50 YAWDEG= YAW*57.296 

RDEG=R"57.296 
RDDEG=RDOT*57 .296 
DDEG=D*57.296 
YAWC=YAWC*57.296 


C WRITE..(6,100)) TIME (XP. XS DOL. vi em 
C lL ,UC,U,UDOT,V,VDOT,YAWC ,YAWDEG ,RDEG , RDDEG , DDEG 
100 FORMAT(IX, TIME= ,F8-3, SEC XP= {hl0.2) Sea 
1, Xe eee 


1' FT XDOT=",F8.4,' PI/SEG) Y=’ F8.2)  Unheeonn 
1 ,F8.4,' FT/SEC 
1',/,2%,' UC=',F8.4,’ FT/SEC U=" (Ree 4. Er see 
iL oor? LO). 
1 ' FT/SEC**2 v=" ,F8.4,' FI/SEC  VDOT="—Rld@ece 
1" FI/SEC=*2' ,/,2X5 YAWG= -BSe4.) DEG seus 
1 ,F15.7,° DEG YAW RATE=". Pl5e7 3) DEG Ee 
LL eeACCEL—« 
1 ,F15.7,' DEG/SEC**2' ,/,2X%, RUDDER] 5 Fl 5 =e 
1 DEG. 
ICOUNT=1 

C TEST IF WANT TO STOP 

300 IF (TIME.GE.ETIME) GO TO 400 

C INTEGRATION STEP SIZE DELT 
DELT=1,0 

C INTEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
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YAW=YAW+R*DELT 
X2=X2+DX2*DELT 
X3=X3+DX3*DELT 
C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
c XDOT=U*DCOS (YAW) -V*DSIN (YAW) 
C YDOT=U*DSIN (YAW) +V*DCOS (YAW) 
C X=X+XDOT*“DELT 
c Y=Y+YDOT*“DELT 
TIME=TIME+DELT 
ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D**2 
GO TO 200 
C J=TDIFF= COST FUNCTION 
400 TDIFF=ISE+ISR 
WRITE(6,500) TDIFF,K1,T1,T2,T3,T4 


500 FORMAT(' ',1X,‘TDIFF =',F15.7,' Kl =',F15.7,' 


Peele PS. 7 52 Xx 
eee. 7 2X, T3= 9, Pi5.7.2xX,'T4=',F15.7) 


REWIND 12 

RETURN 

END 
C BETWEEN LINE 249(END) AND THE FOLLOWING LINE 
C Ve Gono SIN DD -~)WE HAVE TO INCLUDE BOXPLX . 
fm GO.SYSIN DD * 


/ oa 
//GO.FT12FOO1 DD DISP=SHR,DSN=MSS.S2160.A213 


elas 


APPENDIX E 
RESPONSE OF THE SYSTEM FOR IRREGULAR SEAS 


//PROGRA JOB (????,0356),'RESEARCH’ ,CLASS=B 


// EXEC FRTXCLGP, IMSL=DP, REGION=1024K 
//FORT.SYSIN DD * 
C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS 
C HAVE BEEN OBTAINED. 
DIMENSION XX(5) 
C OPTIMAL GAINS FOR CONTROLLER 
XX(1)=2.45967680 
XX(2)=88.2797241 
XX(3)=50.5678864 
XX (4)=5.27039050 
XX(5)=95.3189392 
C THE SUBROUTINE PLANT SIMULATES THE SL-7 CONTAINERSHIP 
CALL PLANT(XX) 
WRITE(6,25) 
25 FORMAT(1X,'OPTIMAL GAINS",/) 
DO 30 I=1,5 
30 WRITE(6,40)1I,XX(I) 
40 FORMAT( 1X, XXC 126) =" - Pala) 
STOP 
END 


C SUBROUTINE PLANT(XX) SIMULATES THE SHIP 
SUBROUTINE PLANT(XX) 
COMMON TDIFF 
REAL*8 1, L2,L3;L4,L5,16 
REAL*8 X,XDOT,Y,YDOT,U,UDOT,V, VDOT, YAW,R,RDOT 
REAL*8 TIME, ETIME,XUDOT, XUU,XVR,XVV,XDD 
REAL*8 YV,YR,YD,YVVR,YVRR, YVVV, YRRR, YDDD, YVDOT 


eZ 


Cl CYS) 


REAL*8 NV,NR,ND,NVVR,NVRR,NVVV,NRRR,NDDD,NRDOT 
Ri omOmeast hy aMZA. ME MASS. DELT 

REAL*8 DYAWE,YAWE,YAWC,ISE,ISR,LAMDA,D 

REAL Sail to. Dawe sDxX2.S5,CHC( LL) ,DX3,xX3,X4 
DIMENSION XX(5) 


CLOSE LOOP ANALYSIS WITH FILTER 


INITIAL CONDITIONS FOR INTEGRATION 
SIMULATION END TIME IN SECONDS 
ETIME=600. 
TIME=0.0 
ICOUNT=1 
INITIALIZE THE COST FUNCTION 
ISE=0.0 
ISR=0.0 
TDIFF=0.0 
civ ere 
GAIN COEFFICIENTS 
K1=XX(1) 
T1=XX(2) 
T2=XX(3) 
T3=XX(4) 
T4=XX(5) 
X,XDOT,Y,YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 


Y=0.0 
XDOT=0.0 
YDOT=0.0 
UeUROT,VY,VDOLT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT=0.0 
VDOT=0.0 
YAW=0.0 
R=0.0 


Heals 


Oec@ 


RDOT=0.0 
YAW=0.0 
ORDERED SPEED IN FEET/SEC 
54.01 FT/SEC=32 KNOTS 
UC=38.81 
AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U=UC | 
D = RUDDER ANGLE 
D=0.0 
L=880.5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*LS5 
SEA DISTURBANCE 
FORCES IN X,Y DIRECTION COMPUTED IN FORCES 
MOMENTS IN Z 
FX=0. 
FY=0. 
MZ=0. 
ISEA IS A SWITCH; ISEA=0(CAL WATER)ISEA=1(SEA STATE) 
ISEA=1 
HYDRODYNAMIC COEFFICIENTS ARE INSERTED HERE AS 
PARAMETERS. 
RHO=1.9876 
MASS=(.0044)*(.5*RHO*L3) 
IZ=(0.00028)*(.5*RHO*LS) 
YAWE=0.0 
X2=0.0 
DX2=0.0 
X3=0.0 
DX3=0.0 
X4=0.0 
200 CONTINUE 
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A C262): ©) 


© 1 @ © 


S=DSQRT(U**2 + V**2) 
INPUT YAW COMMAND 
YAWC=0.0 
IF (TIME.GE.0.0) YAWC=0.0 
ERROR SIGNAL TO DRIVE RUDDER(YAW ACTUAL-YAW ORDERED) 
( COMPENSATOR FILTER ) 
YAWE=YAW - YAWC 
DX2= (YAWE-X2)/T2 
X4=K1*(T1*DX2+X2) 
DX3=(X4-X3)/T4 
D=(T3*DX3+X4) 
AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 


XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 

FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 

XUDOT=(- .0001)*(.5*RHO*L3) 
XUU=(-0.0003)*(.5**RHO*L2 ) 

XVR= (0.0039 )*(.5*RHO*L3) 
XVV=(- .0012)*(.5*RHO*L2 ) 
XDD=(-0.0005 )*(.5*RHO*L2*S**2 ) 

LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 
YV=(-0.00758 )*(.5*RHO*L2*S) 
YR=(0.0023)*(.5*RHO*L3*S ) 
YD=(0.00145 )*(. 5*RHO*L2*S**2 ) 
YVVR=(0.01)%*(.5**RHO*L3/S) 

YVRR=(-0.008 )**(.5**RHO*“L4/S) 
YVVV=(-0.03)*(.5**RHO*L2/S) 
YRRR= (0.003 )*(.5**RHO*L5/S) 
YDDD=(-0.0005)*(.5*RHO*L2*S$**2 ) 
YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 


YVDOT=(-0.0039)*(.5**RHO*L3) 
YVDOT=-2304300.00 
MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 


ele 


qaqa Qa 


30 


35 


QO) 


NV=(-0.00213)*(.5*RHO*L3*S ) 
NR=(-0.00105)*(.5*RHO*L4*S ) 
ND= (-0.0007 )*(.5*RHO*L3*S**2 ) 
NVVR=(-0.015)*(.5*RHO*L4/S) 
NVRR=(-0.008)*(.5*RHO*L5/S) 
NVVV=(0.01)*(.5*RHO*L3/S) 
NRRR=(-0.006)*(.5*RHO*L6/S) 
NDDD= (0.0001)*(.5*RHO*L3*S**2 ) 
NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 


NRDOT=(-0.00027)*(.5*RHO*LS ) 
NRDOT=-1.5096E+11 
SETS SEA STATE TO ZERO 
IF (ISEA.EQ.1) GO TO 30 
FX=0. 
FY=0. 
MZ=0. 
GO TO 35 
UNIT 12 HAS THE SEA STATE DATA NAMED CH 
IT MUST BE SYNCHRONIZED BY TIME 
READ (12) CH 
FX= CH(3) 
FY= CH(4) 
MZ= CH(8) 
CONTINUE 
U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=- XUU*UC**2 
EQUATIONS OF MOTION 
UDOT=( (XVR + MASS)*V*R + XUU*U**2 + XVV*V%%2 
1 + XDD*D*D + FX + XP )/(MASS-XUDOT) | 
VDOT=(YV*V + (YR-MASS*S)*R + YD*D + YVVR*V**2*R 
eeayRR ye Rog 


iG 


L + YVVV*V%"3 + YRRR*R**3 + YDDD*D**3 
1 + FY )/(MAS-YVDOT) 
RDOT=( NV*V + NR*R + ND*D + NVVR*V2*2*R 
lL + NVRR*V*¥R**2 
lL + NVVV*V**3 + NRRR*R**3 + NDDD*D**?*3 
1 + MZ )/(IZ-NRDOT) 
C WHEN TO PRINTOUT 
IF (ICOUNT.EQ.2 ) GO TO 50 
GO TO 300 
C CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW*57.296 
RDEG=R“57 .296 
RDDEG=RDOT*57 .296 
DDEG=D*57.296 
YAWC=YAWC* 57.296 
WRITE (6,100) TIME, YAWDEG 
100 FORMAT(1X,F12.8,1X,F12.8) 
ICOUNT=1 
Ge TeST IF WANT TO STOP 
300 IF (TIME.GE.ETIME) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT=1. 
C INTEGRATION 
U=U+UDOT*DELT 
V=V+VDOT*“DELT 
R=R+RDOT“DELT 
YAW=YAW+R*DELT 
X2=X2+DX2*DELT 
X3=X3+DX3*DELT 
C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT=U*DCOS ( YAW)-V*DSIN (YAW) 
YDOT=U*DSIN(YAW)+V*DCOS (YAW) 
X=X+XDOT*“DELT 
Y=Y+YDOT*“DELT 
TIME=TIME+DELT 


EV7 


ICOUNT=ICOUNT+1 
ISE=ISE + LAMDA*YAWE**2 
ISR=ISR + D**2 
COON Z00 

C J=TDIFF= COST FUNCTION 

400 TDIFF=ISE+ISR 
WRITE(6,500) ISE,ISR,TDIFF 

500 FORMAT('1',5X,'ISE=',F15.7,' ISR=',F15.7,' 
1 TOTAL=',F15.7) 
STOP 
END 

//GO.SYSIN DD * 

/* 

//GO.FT12F001 DD DISP=SHR,DSN=MSS.S2160.A211 
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